JOURNAL OF GEOPHYSICAL RESEARCH: SPACE PHYSICS, VOL. 119, 1-17, doi: 1 0.1 002/20 13JA0 19282, 2014 


Open boundary particle-in-cell simulation 
of dipolarization front propagation 

Alex Klimas, 1 ’ 2 Kyoung-Joo Hwang, 1 Adolfo F. -Vinas, 3 and Melvyn L. Goldstein 3 

Received 2 August 2013; revised 6 November 2013; accepted 9 December 2013. 

[i] First results are presented from an ongoing open boundary 2 'AD particle-in-cell 
simulation study of dipolarization front (DF) propagation in Earth’s magneto tail. At this stage, 
this study is focused on the compression, or pileup, region preceding the DF current sheet. We 
find that the earthward acceleration of the plasma in this region is in general agreement with a 
recent DF force balance model. A gyrophase bunched reflected ion population at the leading 
edge of the pileup region is reflected by a normal electric field in the pileup region itself, rather 
than through an interaction with the current sheet. We discuss plasma wave activity at the 
leading edge of the pileup region that may be driven by gradients, or by reflected ions, or both; 
the mode has not been identified. The waves oscillate near but above the ion cyclotron 
frequency with wavelength several ion inertial lengths. We show that the waves oscillate 
primarily in the perpendicular magnetic field components, do not propagate along the 
background magnetic field, are right handed elliptically (close to circularly) polarized, exist in 
a region of high electron and ion beta, and are stationary in the plasma frame moving 
earthward. We discuss the possibility that the waves are present in plasma sheet data, but have 
not, thus far, been discovered. 
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dipolarization front propagation, J. Geophys. Res. Space Physics, 119, doi:10.1002/2013JA019282. 


1. Introduction 

[ 2 ] It is generally accepted that a dipolarization front (DF) in 
Earth’s magnetotail is a north-south oriented, thin dawn to 
dusk directed current sheet with associated disturbances that 
can propagate large distances, usually earthward, as a recog- 
nizable coherent structure [Rimov et al., 2009; Sergeev et al., 
2009; Rimov et al., 2011a; Sergeev et al., 2012], The current 
sheet, whose thickness remains of the order of a thermal 
proton gyroradius [Hwang et al., 2011; Schmid et al., 201 1], 
separates two magnetic domains with a significantly stronger 
northward directed magnetic field, B z , on its trailing side 
and, often, a small dip in this field component just preceding 
it. Thus, an identifying characteristic of a DF in spacecraft 
observations is the resulting asymmetric bipolar structure in 
the B : tune series with the following increase in the field 
strength appearing almost step-like. In the magnetotail, DFs 
are often observed at the leading edges of fast earthward flows 
and flow bursts within bursty bulk flow events [Angelopoulos 
et al., 1994; Rimov et al., 2009; Hwang et al., 2011; Runov 
et al., 2011a, 2011b]. As such, DFs are thought to separate 
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the cold, dense, ambient plasma sheet population from the 
hot, tenuous outflow of reconnection activity deeper in the tail. 
It has been noted that DFs propagate at approximately the local 
Alfven speed, which is consistent with this scenario [Hwang 
et al., 2011]. Further, particle-in-cell (PIC) simulations of 
impulsive reconnection events have supported this interpreta- 
tion [Sitnov et al., 2009; Sitnov and Swisdak, 2011]. The fast 
flows slow rapidly as they encounter Earth’s dominant dipolar 
magnetic field and a strong pressure gradient [Shiokawa et al., 
1997; Dubyagin etal., 2010; Xing etal., 2010] in the transition 
region between the dipole-like inner and tail-like outer plasma 
sheet. It is thought that substorm or substorm-like auroral 
breakup activity is triggered in the inner transition region by 
this braking mechanism [Ohtani, 2004; Nishimura et al., 
2010; Mende et al., 2011; Sergeev et al., 2012]. Thus, the 
flows provide a link between reconnection activity deeper in 
the tail with consequent substonn and substorm-like auroral 
activity. The flow-burst braking mechanism involved in this 
coupling in the inner magnetosphere is not well understood. 
This mechanism and the role of the DF in plasma acceleration 
and injection into the inner magnetosphere are subjects under 
investigation at present[Nakamura et al., 2009; Ohtani et al., 
2009; Sergeev et al., 2009; Dubyagin et al., 2010; Panov 
et al., 2010a, 2010b; Xing et al., 2010; Hwang et al., 2011; 
Sergeev et al., 2012], 

[ 3 ] With the arrival of a DF at the position of a spacecraft in 
the plasma sheet, it has been found that the appearance of the 
earthward flow does not coincide exactly with that of the 
current sheet [Sergeev et al., 2009; Runov et al., 2011b, 
2011a; Sergeev et al., 2012], The current sheet is often 
embedded within the leading portion of the fast flow with a 
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precursor in front of it. This leading “compression” or “pileup” 
region typically lasts approximately one minute during which 
the flow speed and plasma pressure and density rise slowly 
before the arrival of the current sheet. The pileup region is 
clearly driven forward by the approaching current sheet; it 
appears to consist of a mixture of upstream plasma plus plasma 
that has been reflected forward by the impending DF proper. 
There are two available explanations for the reflection mecha- 
nism: one based on MHD force balance [Li et al., 201 1] and 
the second based on a kinetic model of individual particle 
reflection at the current sheet [Zhou et al., 2010, 2011, 
2012a, 2012b], 

[ 4 ] Li et al. [201 1] have studied a THEMIS [Angelopoulos 
et al., 2008; Sibeck and Angelopoulos , 2008] DF event 
observed in the magnetotail at X — -10 R E by the closely 
separated P3/THD, P4/THE, and P5/THA probes, and at 
X — - 187? £ by the P2/THC probe. Based on the time delay 
between the inner and outer observations, the DF traveled 
earthward at better than 500km/s. The DF was embedded 
within a fast earthward directed bulk flow. Starting approxi- 
mately one minute before the arrival of the current sheet, the 
flow speed slowly increased at the inner probes, consistent with 
the results of Runov et al. [201 lb, 2011a], Using the fortuitous 
placement of the inner probes, Li et al. were able to evaluate 
the magnetic curvature force density and the total pressure gra- 
dient force density acting on the plasma in the neighborhood of 
the current sheet. They found that despite the increased radius 
of curvature following the current sheet, the earthward directed 
curvature force density increased due to the magnetic field in- 
crease within the dipolarized flow burst. Since, in this follow- 
ing region, the pressure gradient force remained unchanged 
compared to quiet conditions, the magnetic field curvature 
force was found responsible for accelerating the plasma in 
the bulk flow region following the current sheet. In the leading 
pileup region, Li et al. found the normally tailward directed 
pressure gradient force reduced compared to quiet conditions 
while the curvature force density remained consistent with 
quiet conditions. Thus, they concluded that the pressure gradi- 
ent force was responsible for the formation of the leading 
pileup region and for the earthward acceleration of its ambient 
plasma. Li et al. suggested that flow bursts accompanied by a 
large amplitude DF will penetrate deeper into the dipole field 
dominated near-Earth plasma sheet, implicating the DF in the 
flow burst braking and auroral breakup triggering process. 
Zhou et al. [2010, 2011] had earlier shown that the formation 
of the compression region could be explained as an ion pickup 
region populated by particles reflected at the current sheet. Li 
et al. suggested that a distribution function analysis could tell 
which of the two mechanisms was responsible for the appear- 
ance of the region. 

[ 5 ] In a series of recent papers, Zhou et al. [2010, 2011, 
2012a, 2012b] have investigated the relationship between 
fast flows, DFs, associated pileup regions, and the emergence 
of field aligned particle beams in the plasma sheet boundary 
layer (PSBL). Their results follow from their case studies and 
statistical surveys of THEMIS 2008 and 2009 tail season 
observations plus supporting test particle simulations. In its 
latest version, the simulation model is based on a realistic 
two-dimensional kinetic plasma sheet model, a slight modifica- 
tion of the Pritchett and Buchner [1995] model, containing a 
monotonically decreasing neutral sheet B : with increasing depth 
into the tail. On top of this background, an electromagnetic 


model of a propagating DF is imposed. This consists of a 
B : field with a finite width in the dawn-dusk direction and 
containing a sharp step-like increase at a position that prop- 
agates earthward at a fixed speed. A dawn-dusk inductive 
electric field is also included to satisfy Faraday’s Law for 
the varying magnetic field. The propagating DF model has 
no particle content. At each step in the simulation, taking 
the presence of the DF component into account, particles 
are integrated backward in time from their position at a vir- 
tual spacecraft to their origin in the background plasma 
sheet model. Liouville’s theorem is invoked to relate the 
particle distribution at the spacecraft to that at the origin 
in the kinetic background model. The result is the forma- 
tion of an energetic particle precursor to an arriving current 
sheet that is populated by particles that have been reflected 
and accelerated at the current sheet. The model does repro- 
duce observed upstream earthward and dawnward drifts 
[Zhou et al., 2012b], Further, in this model, particles in 
the central plasma sheet arrive 10’s of seconds early while 
particles in the PSBL arrive minutes early, in very good 
agreement with THEMIS observations [Zhou et al., 
2012a], From this point of view, the pileup region is populated 
by particles reflected and energized by the magnetic and 
electric fields of the traveling current sheet. The particles 
in the central plasma sheet are limited to a region earthward 
of the current sheet of extent approximately one thermal 
gyro radius in the B z field while particles in the PSBL move 
more freely earthward. 

[6] In this paper, we report first results from an ongoing 
effort to develop a particle-in-cell (PIC) simulation model 
that will be analogous to the Zhou et al. [2012a] model, but 
fully self-consistent. Our intent is to use newly available open 
boundary techniques [Klimas et al., 2008, 2010, 2012] to 
drive a DF across a simulation domain that initially repre- 
sents a quiet plasma sheet. Initial results that we present 
below demonstrate, first, that the open boundary, driven 
simulation approach is feasible. Further, we find that many 
of the observed DF properties are reproduced by this model. 
A pileup region is produced whose properties are generally 
consistent with the force balance picture [Li et al., 2011] but 
differ from the particle pick up picture [Zhou et al., 2010, 
2011, 2012a, 2012b]. The simulated pileup region is, in- 
deed, populated by gyrophase bunched reflected ions. 
However, these ions are reflected through an interaction 
with the pileup region itself, not through an interaction with 
the current sheet. Such a self-consistent interaction possi- 
bility is not included in the Zhou et al. particle tracing 
model. The present version of our model is highly ideal- 
ized. It will be important to verify these preliminary results 
using more sophisticated versions of the driven DF PIC 
simulation approach. 

[ 7 ] We find that wave activity at the leading edge of the 
pileup region is a persistent feature of our DF simulations. 
Wave activity growth coincides with the maximum gradient 
in density, pressure, field strength, etc., at the leading edge. 
We show that this is also the region where the ion velocity 
distribution is highly nongyrotropic due to the presence of 
reflected ions. Both the steep gradients and the nongyrotropic 
velocity distributions are candidates for exciting the waves. 
However, we have not identified the wave mode(s) involved. 
The possibility that the waves are present in plasma sheet 
data but have not, thus far, been discovered is discussed. 
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Figure 1. Driven DF simulation setup. 


2. Simulation Setup 

[s] Figure 1 gives an overview of the simulation setup for 
the driven DF model. The simulation models the passage of 
a propagating DF through a portion of the central plasma sheet 
in Earth’s magnetotail. The earthward direction is to the left, in 
the direction of increasing x, with the deeper tail to the right. 
The z axis points to the north and the y axis, from dawn to 
dusk. Plasma is driven into the simulation domain through 
the open right-side boundary through the application of an 
out-of-plane electric field, E v , on the boundary and plasma is 
free to pass out of the simulation domain through the open 
left-side boundary. The initial magnetic field is solely in the 
z direction. In this version of the model, the upper and lower 
boundary conditions are periodic. This, plus the lack of an x 
component in the initial magnetic field, produces the most 
serious limitations in the applicability of this idealized model. 
Work on a more generalized model in which these restrictions 
are removed is in progress. It will be very important to verily 
that the conclusions that we present in this paper remain valid. 

[ 9 ] The simulation is initialized with the magnetic field in 
the z direction (northward) given by [ Pritchett , 2008], 

B, = B 0 (r- tanhQ) (1) 

with a current sheet in the y direction (dawn-dusk) of thick- 
ness A in the x direction given by 

^ = T cosh ~ 2 © (2) 

[ 10 ] The parameter B 0 measures the initial overall strength 
of the DF, A can be adjusted to control the initial thickness of 
the DF, and R is a free parameter that can be used to control 
the relative strengths of the field on either side of the current 
sheet. Notice that 

B : ~ 5 0 (fl+l) 

x/ (3) 

B z ~ B 0 {R- 1) 

x/AM) 

[ 11 ] Standard symmetric reconnection simulations use 
R = 0, asymmetric simulations (e.g., Pritchett [2008]) are 


based on 0 < R < 1 , and for this case we choose R > 1 , with 
the initial ratio of the field strengths on either side of the 
DF controlled by our choice. 

[ 12 ] To achieve an equilibrium initial state, we balance the 
magnetic field pressure with the plasma pressure so thatP(x) + 
S?/8w = const = C. We assume Gaussian particle distribu- 
tions in this initial state and further assume that P(x) = n(x)k 
( Tj+T e ), in which n(x) is the ion or electron density and the 
ion and electron temperatures are unifonn in x. We make the 
further assumption that B\/%k = anok(T, + T e ) in which a 
is a free parameter and n o = n(x=0). Finally, then, we have the 
particle density 

n(x ) = «o(l — 2at?tanh 

which can be used to load the particles in the simulation code 
to create the initial state. Notice that 

»(*) «o[l -«(1 -2fl)] 

1 (5) 

n(x) ~ w 0 1 — «(1 +2tf)l 

x/X>0 1 1 

and therefore 


to keep the density nonnegative on the low density side. This 
last condition limits the energy density of the field relative to 
that of the particles. Since we must have R > 1 to avoid a null 
in the magnetic field, a< 1/3. 

[ 13 ] All variables in the DF simulation code are dimension- 
less. To relate their values to a dimensional environment, two 
fundamental parameters must be specified, n 0 and B 0 . Both of 
these relate to the initial state of the simulation. Assuming no 
dependence on the z coordinate, the particle density (equal 
for ions and electrons) at the central position in x within the 
current sheet is given by n 0 and the magnetic field strength at 
the same position and time is RB 0 . In applications to observed 
DF data, the field strength and particle density should be set to 
values observed in the central plasma sheet portion of the DF 
current sheet. Using the symbol * to denote dimensionless 
simulation variables, then the dimensional density n = n 0 n 
and, similarly, the magnetic field B=i> 0 B , time t =a> ic t in 
which w ic is the ion cyclotron frequency in the B 0 field 
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- 10 




X Coordinate (arbitrary origin) 

Figure 2. B z at indicated times (blue). (Curve at / = 1 9 has 
been eliminated for clarity. B- strength in back of current 
sheet does not vary significantly from 7= 15 to 7= 19.) 
Dawn-dusk current density at indicated times (red). All 
curves show north-south averages over z. 

strength, lengths are normalized with respect to the ion inertial 
length d i =c/o)p i =c(4jie 1 no/m l )~ 12 (cgs units), particle and 
moment velocities are expressed as v = V A0 \ , in which V A0 
is the Alfven speed computed using the ion mass and n 0 and 
B 0 for the density and field strength, and the electric field 
E = V A oB 0 E . The dimensional and dimensionless ion and 
electron pressure tensors are related through P = (B^/4n)V* 
and the temperatures through T* = (An / BfynqkT . 

[ 14 ] We have used the 2 lA-dimensional electromagnetic 

PIC code described in Hesse et al. [1999], modified to incor- 
porate the open boundary conditions [Klimas et al., 2008, 
2010, 2012], The electromagnetic fields are integrated implic- 
itly to avoid the Courant constraint on the propagation of light 
waves. The simulation discussed here (denoted DF-1 in the 
following) was initialized with approximately 4 x 1 0 8 particles 
on a 1025 x 625 grid in the x x z directions. The system size is 
L x = 65.6dj by L z = 40c/,. Two particle species, one of ions and 
one of electrons, with mass ratio mjm e = 25 were included. As 
a matter of simplicity in this first-of-its-kind simulation, 
particle temperatures were initialized with T t =T e . This 
choice is in keeping with the intent to drive the DF long 
enough to induce a loss of memory of the initial state of 
the simulated plasma. We show below, however, that the 
simulated DF propagates at super Alfvenic speed and, 
consequently, interacts with upstream plasma that remains |> 
unchanged in the initial state. In future studies, we will initial- > 
ize with the more typical magnetotail 7) >57),. It will be impor- J 
tant to explore the significance of this change. Using the ^ 
electron cyclotron frequency co ec = eBo/rngC and the electron § 
plasma frequency co ep = (Ane 2 n 0 lm e ) 112 , we set to ep lw ec =A. Jj 
To ensure particle energy conservation, very small integration “ 
time steps were used; At =0.0005 in the dimensionless time g 

defined in the preceding paragraph. The initial current sheet d 

thickness was set to 0.5(7,, the parameter 7? = 1.5, and a = 0.2. 

The simulation was driven through the application of an out- 
of-plane electric field component E y at the inflow boundary. 

This was ramped up from its initial value, 0, to its final value, 

15, over the time interval, 7=0 — 5. 

[ 15 ] Charge conservation is guaranteed by an iterative appli- 
cation of a Langdon-Marder type correction [Langdon, 1992] 


to the electric field at each time step. At the open boundaries, 
we apply a gradient free boundary condition on the out-of-plane 
electric field component, set the normal component to zero, 
^ and we leave the z-component unchanged from the previous 
o> time step. Then, the Langdon-Marder correction is applied to 
g correct the electric field over the entire grid, including at the 
% boundaries. For the magnetic field at the open boundaries, 
o we set the out-of-plane magnetic field component to zero. 
£ At each step in the iterative implicit field solver, we set 
z the z-component of the magnetic field to satisfy V x B=j 
and then the normal component to satisfy V • B = 0. This 
totality of particle and field boundary conditions generally 
leads to stable, quiet solutions with no trace of charge 
buildup near the boundaries or anywhere else on the 
computational grid. 

3. Simulation Overview 

3.1. Dawn-Dusk Current Density and B. Strength 

[ 16 ] Figure 2 provides an overview of the magnetic field 
and cross-tail current density encompassing most of the 
DF-1 simulation, which evolved from 7 = 0 to 7 = 20. For this 
figure and for the remainder of this paper, quantities denoted 
as north-south averaged will have been averaged over the full 
range of the z-coordinate. It can be seen that both the strength 
of the field jump at the current sheet and the strength and 
width of the current sheet evolve considerably under the 
influence of the driver field. Notice the development of a 
precursor, a pileup region, in the field strength. The simula- 
tion was stopped as the current sheet approached the outflow 
boundary and the majority of the pileup region had already 
exited. The continuing growth of the current sheet intensity, 
as well as the appearance of a bifurcation of the current sheet 
(inset), shows that the DF never came into equilibrium with 
the driver. At 7=0, the current sheet was stationary and then 
only accelerated in response to the driver field later. As it 
accelerated, a pileup region began to form and then grow. 
Late in the simulation, the pileup region extended too far in 
front of the current sheet to be a valid model of the phenom- 
enon in the magnetotail. This growth is a direct consequence 
of the combination of periodic boundary conditions in the 
north-south directions plus the lack of an x component in 



Time 


Figure 3. Plasma inflow speed at inflow boundary (green). 
Propagation speed of current sheet peak (red). Local Alfven 
speed at the position of the current sheet peak (blue). 
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Figure 4. Showing the development of the pileup region at 
t=l and t= 15. Vertical dashed lines show the positions of 
the current sheet at these times. Ion pressure tensor compo- 
nent P xv (blue). Ion earthward flow speed (red). All quantities 
north-south averaged. 


3.3. Pressure and Flow Speed in Pileup Region 

[is] The compression, or pileup, region preceding the 
arrival of a DF current sheet is typically detected in the growth 
of the plasma pressure and flow speed before the arrival of the 
current sheet [Runov et al., 2011b, 2011a; Sergeev et al., 
2012]. Figure 4 shows snapshots at two times of the ion 
pressure tensor component P xx and the ion flow speed Re- 
taken from the DF-1 simulation. The positions of the current 
sheet at these two times are indicated by the vertical dashed 
lines. It can be seen that the spatial extents of both the pres- 
sure and flow speed precursors increase with increasing 
time while the pressure (and density, not shown) increases 
in amplitude as well. The relatively low pressure in back 
of the DF is due to the low particle density in that region. 
At the later time, the DF and local plasma in front and in 
back of the DF all move together at the same speed, in 
agreement with recent THEMIS spacecraft observations 
[Runov et al ., 2009; Sergeev et al., 2009; Zhou et al., 
2011], At the earlier time, this is not exactly true but it can 
be seen from Figure 3 that the DF was still being accelerated 
somewhat at that time. 


the initial magnetic field configuration. Plasma remained 
trapped in front of the propagating current sheet with no 
escape route, as in the magnetotail. In the next section, we 
show that by t = 7 most of the acceleration of the current sheet 
had finished. From then until approximately t= 13, the simu- 
lated pileup region remained an acceptable model of the mag- 
netospheric counterpart. We have used simulation output 
from before and after those times in limited ways discussed 
below but our conclusions concerning the properties of the 
simulated pileup region are limited to its behavior in the 
interval 1 <t< 13. 

3.2. Flow and Alfven Speeds 

[ 17 ] Plasma was driven into the simulation domain through 
the imposition of a uniform out-of-plane electric field at the 
boundary. The inflow speed is given by EJB- at the boundary. 
The green curve in Figure 3 shows this inflow speed. The initial 
growth in speed up to 7=5 is a consequence of the ramping up 
of E y during that interval. The following slow decline, with 
asymptote somewhat above 2, reflects the growth with increasing 
time in B z at the boundary, which is evident in Figure 2. The red 
curve shows the speed at which the peak of the current sheet 
traveled toward the outflow boundary. After some delay while 
the driver field propagated to its position, the current sheet 
accelerated rapidly. The decline in this speed near the end of 
the simulation would be temporary in a longer simulation; the 
red and green curves must eventually come into agreement. 
The blue curve shows the evolution of the local Alfven speed 
at the position of the current sheet peak. The large jumps in this 
local speed are due to the combination of a very steep gradient 
in B- strength at the current sheet position and a discrete grid on 
which to measure the strength. The local Alfven speed appears 
to grow approximately linearly with time, although the noise 
near the end of the simulation admits the possibility that the 
speed may be saturating, as might be expected from the appar- 
ent saturation of the B- strength in back of the current sheet 
evident in Figure 2. With this proviso, we find that both 
current sheet speed and the local Alfven speed self-adjust to 
approximate equality, in agreement with the observations of 
Hwang et al. [2011] in the central plasma sheet. 


3.4. Current Sheet Details 

[ 19 ] Observed DF current sheet thicknesses are typically of 
the order of 1-3 thermal ion gyroradius or ion inertial lengths 
(see Sergeev et al. [2012] and references therein). Panel (a) of 
Figure 5 provides a comparison of these length scales, at 



Figure 5. (a) Current sheet density (red). Thermal proton 
gyroradius in current sheet neighborhood (blue), (b) Charge 
density (red). Normal electric field component E x (blue). 
dEJdx (green). 
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Figure6. Nonideal E*=E X + [V e x B| v (blue). Hall contribu- 
tion to E* from (8) (red). Summed Hali and pressure gradient 
contributions from (8) (green). 


t= 10, for the DF-1 simulation. The x coordinate is measured 
in units of the proton inertial length, computed using the 
proton density in the current sheet at t= 0. Since the proton 
density in the current sheet does not change significantly 
over the entire simulation, the x coordinate is effectively 
measured in units of the local ion inertial length at this time. 
The red curve in panel (a) gives a detailed view of the cur- 
rent sheet density. Depending on how it is measured, it 
can be seen that the current sheet thickness is of the order 
of a proton inertial length, or less. The blue curve shows 
the thermal proton gyroradius in the neighborhood of the 
current sheet. Because of the jump in field strength at the 
position of the current sheet, the gyro radius varies from 
— 3.5 in front of the current sheet to — 0.3 in back. At the 
position of the current sheet peak, the gyro radius is — 0.7. 
Thus, we find that the inertial length, gyroradius, and cur- 
rent sheet thickness are all approximately equal. 

[ 20 ] The steep growth in magnetic field strength at the cur- 
rent sheet position acts as a barrier that prevents the denser 
upstream plasma from entering the downstream region. 
Because the upstream proton gyro radii are comparable to 
the thickness of the current sheet, they are able to penetrate 
the magnetic barrier at the current sheet to a depth that is 
comparable to the current sheet thickness while the electrons 
are turned back more quickly. The consequent separation of 
electrons and protons leads to a layer of charge coincident 
with the position of the current sheet. The red curve in panel 
(b) of Figure 5 shows this charge distribution at t= 10 in the 
DF-1 simulation. The negative peak to the left gives an indi- 
cation of the penetration depth of the electrons while the 
positive peak to the right indicates the proton penetration 
depth. Consistent with the charge distribution, an earthward 
directed electric field E x is produced in a thin layer at the 
current sheet position. The blue curve in panel (b) illustrates 
this field. The positive pulse joins the positive and negative 
charge layers and is effectively shielded from the upstream 
region. The positive charge layer is unshielded to the right, 
leading to a tailward directed E x in the downstream region. 
The green curve in panel (b) is the result dEJdx, confirming 
that the field is the direct result of the separation of electrons 
and protons at the magnetic barrier. Consequently, E x might 
be considered a Hall type field. 


[ 21 ] DF observations often exhibit a strong earthward 
directed Hall type electric field normal to and localized at 
the current sheet [Zhou et ah, 2009; Runov et al., 2011a; 
Sergeev et al, 2012], Zhou et al. [2009] have discussed a 
THEMIS event that exhibited a strongly peaked earthward 
directed DC electric field normal to the DF in the plasma rest 
frame, i.e., a DC peak in the nonideal E* = E x + (V x B) x 
component. They found that the peak could be very well fit 
by the Hall term (j x B )Jne and, thus, concluded that the 
DC field was a Hall field. Runov et al. [201 la], however, in 
a superposed epoch study have found that the median E* is 
well explained by the Hall term in front of the DF but not 
in back. For the DF simulation under discussion here, we find 
that the electric field of Figure 5b has two equally important 
contributions to its nonideal part, the Hall contribution and an 
opposing electron pressure gradient contribution. 

[ 22 ] The nonideal corrections to the V x B drift contri- 
bution to E x can be obtained from the x-component of 
the dimensionless electron or ion momentum equations. 
Neglecting other than the ^component of the magnetic 
field, these can be written 


E* — E x + [V„ x B] X —E x + V ay B- 

1 

n a q , 


§ t (m a n a r ae )+l{m a « a V„r m + P ] „) 


(7) 


in which E* is the nonideal portion of the electric field, the 
symbol a stands for either electron or ion quantities with 
q a = ± 1, and P a is the electron or ion pressure tensor. For 
electrons in the DF simulation, we find that the thermal pres- 
sure gradient term on the right side of (7) dominates strongly. 
Thus, assuming Ui = n e , we can write 



1 8 

n e 8x 


( 8 ) 


[ 23 ] The blue curve of Figure 6 shows the nonideal compo- 
nent of the E x illustrated in Figure 5. The contribution of the 
Hall term on the right side of (8) to E* is shown by the red 
curve; it peaks at approximately the value 50. The summed 
contributions of both terms on the right side of (8) are given 



Figure 7. Ion V y (red). Electron V v (green). E x B drift 
velocity (blue). Electron V y with pressure gradient correc- 
tion (magenta). 
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Figure 8. Left column: Ion dawn-dusk current distributions at t=l. Right column: Reduced ion perpen- 
dicular velocity distributions for particles collected from regions indicated by green rectangles superposed 
on current distributions to their left. 


by the green curve. Since it is precisely the difference between 
n, and n e that is responsible for generating E x , we have con- 
firmed that this result is insensitive to the choice of «, or n e 
in the denominator of the Hall term, or even to the choice 
( jy/ n i) + ( j y Jn e ) =V iy — V ey . Thus, we conclude that in this DF 
simulation the nonideal electric field normal to the current sheet 
contains a large Hall contribution but also an almost equally 
large, and opposing, electron pressure gradient contribution. 

3.5. Dawnward ion and Electron Fluxes in 
Pileup Region 

[ 24 ] In a pileup region, near the magnetotail neutral sheet, 
ion fluxes are generally observed in the dawnward direction 
[Hwang et al., 2011; Zhou et al., 2011, 2012b]. Zhou et al. 
[2012b] have explained this behavior using the theory of 
Buchner and Zelenyi [1989] for particle trajectories in the 
magnetotail. From this theory, they were able to show that, 
after being reflected by the DF, only ions that quickly turn 
dawnward with reduced earthward speed can remain near 



Figure 9. Two reduced ion velocity distributions just in 
front of the pileup region at (a) t = 5 and (b) t = 15. 
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Figure 10. At selected times, v 2 values at peaks in reflected ion velocity distributions (red). Least squares 
linear fit (blue). 



the neutral sheet; other ions that retain high earthward speed 
must move toward the PSBL. The result is a preponderance 
of dawnward moving ions near the neutral sheet in the 
pileup region. 

[ 25 ] In the results of the DF-1 simulation, we find a similar 
dawnward ion drift in the pileup region. Below, we show that 
this drift is the expression in the velocity moment of the 
motion of a group of gyrophase bunched reflected ions as they 
begin to gyrate about the B z field after interacting with the DF. 
This behavior appears in the ion velocity moment component 
V y in the neighborhood of the DF, as shown by the red curve of 
Figure 7 at t = 7. The negative excursion, indicating dawnward 
motion in the pileup region on the earthward side of the current 
sheet is evident. Since, as we show below, the electrons do not 
participate in the gyrophase bunched motion of the ions, it is 


somewhat surprising that the electron cross-tail drift in the 
pileup region, as shown by the green curve of Figure 7, is 
almost identical to that of the ions. The small difference that 
can be seen in the figure will be shown due to wave activity 
in the pileup region that involves local current. The blue curve 
of Figure 7 shows the cross-tail E x B drift speed. It can be 
seen that, despite the fact that the ions and electrons are 
drifting at essentially the same speeds, they do not drift in 
agreement with the E x B drift. For the electrons, this disagree- 
ment can be understood as due to the important thermal pres- 
sure gradient term discussed in the preceding paragraph. 
Retaining this term on the right side of (7), we have 



Figure 11. Selected portions of the north-south averaged ratios at t= 15 (red), t= 10 (green), and t= 5 (blue). 
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Figure 12. (a) Time derivative of ion momentum as a func- 
tion of position x (black). Total force density obtained 
through evaluation of all terms in the right-hand side of 
(13) (red). Electromagnetic force density (green). Pressure 
gradient force density (blue), (b) Time derivative of ion mo- 
mentum as a function of position x (black). Electromagnetic 
force density (green). Thermal pressure gradient force 
density (red). Dynamic pressure gradient force density (blue). 
All at 1=7. 

[ 26 ] The magenta curve of Figure 7 shows an evaluation 
of the left side of (9), in good agreement with the E x B drift, 
i.e., the right side of (9). We have found similar results for the 
ion drifts. However, in that case, all of the terms on the right 
side of (7) play important roles and none can be neglected. 

[ 27 ] We are unaware of any published observations of 
dawnward electron flux in the pileup region. From Figure 4 
of Hwang et al. [201 1], it can be seen that the ion dawnward 
drifts are ~ 200 km/s. If the electron drift equaled this value, 
then it would have to be detected in a velocity distribution 
whose thennal width is typically several hundred times as 
large, a very difficult task. We suspect that this difficulty is 
the reason for the lack of published dawnward electron 
drift observations. 

4. Reflected Ions 

[ 28 ] Zhou et al. [2010, 2011, 2012a, 2012b] have shown 
that many observed properties of the DF pileup region can 
be explained using their particle tracing model, which invokes 
ion reflection and acceleration at the DF current sheet. The 
DF-1 simulation exhibits ion reflection and acceleration as 
well; however, as we show in this section, the interaction 


in this case is between ions and the pileup region itself, a 
possibility that is absent from the particle tracing model. 

[ 29 ] The second column of Figure 8 contains examples of 
reduced (integrated over vn = v z ) ion perpendicular velocity 
distributions obtained from the DF-1 simulation at 7=7. The 
images of the first column are all identical; they show the 
dawn-dusk ion current distribution at this time. The dawnward 
ion drift discussed above produces the region of negative ion 
current indicated by the red band to the left of the primary cur- 
rent sheet. Superposed on each image is a green rectangle 
showing the region in the ion current distribution from which 
the ions of the corresponding velocity distribution to the right 
were obtained. All of the distributions are effectively averaged 
over z. The first distribution (b) was constructed using ions 
from a region just in front of the pileup region. It can be seen 
that it consists of an essentially gyrotropic core that is slightly 
displaced from zero drift, plus a gyrophase bunched group of 
energetic, reflected ions in the ( v x , — v y ) quadrant. The second 
distribution (d) shows that this group of reflected ions is 
strengthened deeper in the pileup region, producing, despite 
the shift of the distribution core toward the + v v direction, the 
net dawnward current. Finally, the third distribution (f), from 
the DF current sheet, is roughly gyrotropic, heated relative to 
the distribution in front of the pileup region, and significantly 
shifted in the (v x ,v v ) direction to accommodate the earthward 
propagation of the duskward directed current sheet. 

[ 30 ] The energy of the reflected ions in front of the pileup 
region can be estimated, thereby showing that the particles 
are reflected by a normal electric field in the earthward direc- 
tion and not through an interaction with the current sheet. 
Figure 9 shows two reduced ion velocity distributions from 
positions just in front of the pileup region at two times, the 
first quite early in the simulation and the second quite late. 
The upward shift in speed of the reflected ions with increasing 
time is evident. We have estimated the speed of the reflected 
ions by finding the peak in the distribution of reflected ions. 
Estimates of v 2 obtained in this manner are shown in 
Figure 10 for selected times. The dashed blue line is a 
least-squares fit given by 9.56 + 2.19 x t. 
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X Coordinate (arbitrary origin) 


Figure 13. Dawn-dusk current distribution with north- 
south averaged portion removed. Dashed white line marks 
position of maximum current filament wave activity. 
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Figure 14. North-south averaged current filament wave 
magnitude (red). North-south averaged B z on expanded 
scale (blue). 


[ 31 ] All of the points shown in Figure 10 were obtained 
from velocity distributions just upstream of the pileup region 
that contained, as in Figure 9, undisturbed core distributions 
at the same temperature. The increase in reflected ion v 2 shown 
in Figure 1 0 is a real increase; it does not reflect a change in the 
starting point of the reflected particles. Zhou et al. [20 1 0] have 
shown that reflection of ions at the current sheet should 
increase the speed of the ions by approximately twice the 
speed of the moving current sheet. The approximately linear 
growth with time shown in Figure 10 is incompatible with 
the prediction of this mechanism given the evolution of the 
current sheet propagation speed shown in Figure 3. There, it 
can be seen that the propagation speed has essentially saturated 
by t= 7. If this were the reflection mechanism at work, then the 
upward trend in Figure 10 should also cease at this time. 

[ 32 ] Runov et al. [2011a] have noted, “. . . the dipolarization 
front may be described as a magnetopause of the energetic 
plasma intrusion into the quiet plasma sheet.” In the following, 
we present evidence that suggests the leading edge of the 
pileup region behaves in some ways as a normal shock 
leading this magnetopause. In particular, the mechanism 
for reflecting and accelerating the gyrophase bunched ions 
at the leading edge may be similar to that which operates 
at Earth’s bow shock to produce the upstream gyrophase 
bunched ions observed there [Eastman et al., 1981; 
Gurgiolo et al., 1981, 1983]. 

[ 33 ] In terms of the dimensionless variables in use here, 
and using north-south averaged electromagnetic field compo- 
nents, the motion of an ion upstream of the current sheet is 
governed by 


dv x E x 


d\’y 

dz 


-(v, - V x ) 


(10) 


with z = B z t and in which we have used E y = V X B Z with V x 
being the z-averaged ion velocity moment component. With 
v 2 = + v 2 and treating the field components and V x as 

local constants, (10) yields 


[ 34 ] Leroy etal. [1981] have discussed the results of a 1 -D 
hybrid simulation of Earth’s bow shock in a perpendicular 
configuration. The foot region of their simulated shock is 
analogous to the leading edge of the pileup region in our 
DF-1 simulation. Their explanation for the reflected ions that 
they observed in the simulated shock foot region can be 
stated in terms of equation (10): They found that the electric 
field, E x (expressed as the gradient of a potential), was too 
weak to reflect the incoming ions. Ignoring the electric field 
in (10), we are left with ions gyrating in the B z field about 
the position V x in the velocity space. In their simulation, an 
ion that traveled from upstream into the foot region was 
turned around through this gyration motion while the electric 
field provided the extra acceleration necessary to convert it 
into a reflected ion traveling upstream. In the following, 
we explore the possibility that the reflected ions in the DF-1 
simulation are produced in a similar manner. 

[ 35 ] Consider a group of ions just upstream of the pileup 
region with the highest available v x (v y — 0) in the undisturbed 
thermal velocity distribution. From (10), it can be seen that as 
the pileup region engulfs them these ions will be accelerated in 
the negative y direction and are therefore candidates for the 
gyrophase bunched dawnward drifting ions discussed in the 
preceding section. Figure 11 shows the ratio EJB- in the 
pileup region at three times in the DF-1 simulation. It can 
be seen that the ratio exhibits peaks at the leading edge of 
the pileup region. We will assume that, as in the Leroy 
et al. [1981] simulation, the group of high v x ions are 
reflected through their gyro motion in the pileup region B z 
and accelerated through their interactions with these peaks. 
We will assume that the interactions last for as long as it 
takes to turn the ions around, approximately a half gyro 
period. Then we will find that the energy gained through these 
interactions is consistent with the strengths of these peaks. 




1 c/v 2 

2 dz 


= y Vx + VxVy 


( 1 1 j Figure 15. Evolution of dawn-dusk current density at 
position of maximum current filament wave activity. 
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Figure 16. Current filament wave power distributions on cuts through the peak power position for k z >0 
in terms of (a) inverse wavelength and (b) frequency in units of co ci . 


[36] From (11), a rough estimate of the change in v 2 that 
occurs during one reflection interaction with EJB Z can be 
obtained assuming that the interaction lasts a half gyro period, 
Az = n. Then 

v 2 = J(inii) + 2n(^\vx (12) 

in which W is a measure, to be determined, of this particle ve- 
locity component during the interaction. 

[ 37 ] Figure 1 1 shows the distribution of E X IB : at 1 = 5 (blue), 
7= 10 (green), and 7= 15 (red). To avoid confusion, only the 
portions of the distributions on the leading side of the current 
sheet have been plotted. The green curve, for example, shows 
the extension (divided by B-) into the pileup region of the 
negative E x on the upstream side of the E x peak shown in panel 
(b) of Figure 5. Each EJB- distribution is effectively zero in 
front of a leading peak that coincides (not shown) with the 
leading edge of the pileup region. The increase of the strength 
of this peak with increasing time is evident. If this peak is 
responsible for the reflection of ions at the leading edge of 
the pileup region, then its strength must be compatible with 
the measured growth of v 2 with time shown in Figure 10. 
We use the intercept at 7=0 of the linear fit to the points on 
Figure 10 for v 2 ( init ); i.e., we are assuming that in the limit 
EJB Z — ► 0, v 2 remains at the value v^( init ) and that value is 


given by the intercept value, -10. Further, we use v(init) for 
yT These choices are consistent with the particle lying initially 
on the outer fringes of the core velocity distributions shown, 
for example, in Figure 9. Since, from (11), maximum acceler- 
ation is given to the highest v x available, we should expect to 
find the origin of the gyrophase bunched accelerated ions at 
the highest speed fringe of the core velocity distribution. 
With these choices, we use (12) to solve for the strength of 
EJB- that a particle must interact with for a half gyro period 
to attain the v 2 value given by the linear fit of Figure 10 at 
the times, 5, 10, and 15. Those strengths are indicated by the 
positions of the black crosses on Figure 1 1. It appears that this 
is a plausible explanation for the reflection and acceleration of 
ions at the leading edge of the pileup region; it does not 
involve reflection at the current sheet, as in the Zhou et al. 
[2010, 2011, 2012a, 2012b] model. Notice, from Figure 10, 
the energy of the reflected ions is increased by approximately 
a factor of 2.5 at 7=7. 

5. Force Balance 

[38] Li et al. [201 1], using an appropriate configuration of 
three of the THEMIS spacecraft, have studied changes in the 
total plasma pressure gradient and magnetic curvature force 
densities coincident with the arrival of a DF in the plasma 
sheet. Using the spacecraft configuration to model the 
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Figure 17. Evolution of (a) B x and (b) B v at position of maximum current filament wave activity. 
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Figure 18. Time averages of the magnitudes of the 
magnetic field component perturbations shown in Figure 17. 
Time averaged B x wave amplitude (red). Time averaged B v 
wave amplitude (blue). 

curvature of the magnetic field and assuming a constant prop- 
agation speed to determine spatial gradient scales, they were 
able to quite accurately predict changes in the plasma flow 
speed in the DF neighborhood. In back of the DF current sheet, 
they found that, while the curvature of the field lines de- 
creased, the large increase in field strength led to an increase 
in curvature force that dominated the effects of the pressure 
gradient force. This aspect of their observations is modeled 
in the DF-1 simulation by the large jump in B z at the posi- 
tion of the current sheet. Li et al. found that the plasma in the 
pileup region was accelerated earthward by a reduction in the 
normally tailward directed pressure gradient force while the cur- 
vature force in this region remained relatively fixed. The 
undisturbed plasma in the DF-1 simulation does not contain a 
tailward pressure gradient force or a magnetic curvature force 
to balance it. Nevertheless, we do find that the disturbed pileup 
region plasma is accelerated earthward due primarily to an 
earthward directed pressure gradient force in that region. 

[ 39 ] The x-component of the ion momentum p ix = 
m i n iVi X obeys 

+ (13) 

[ 40 ] There are two driving forces for changing the ion mo- 
mentum. We will denote these as the electromagnetic force 
density (first term on the right side) and the pressure gradient 
force density (second term). The pressure gradient force con- 
tains two sources, the dynamic pressure (bulk velocity com- 
ponents) and the thermal ion pressure (P XXI ). These can be 
evaluated directly from the DF-1 simulation. 

[ 41 ] The black curve in panel (a) of Figure 12 shows a 
direct evaluation, at f=7, of the time derivative on the left 
side of (13). The red curve shows the agreement of a super- 
position of all of the contributions on the right side of (13). 
The peaks in this figure are at the position of the current 
sheet. In back of the current sheet, the approximately uniform 
plasma is moving earthward, leading to no significant change 
in the ion momentum. At the position of the current sheet, 
there is a net negative acceleration of the ions. The green 


curve shows the expected positive acceleration due to the elec- 
tromagnetic (primarily magnetic) force, which is, however, 
overwhelmed by the negative peak in the pressure gradient 
force (blue curve). The negative peak can be understood 
through an examination of Figure 4. There, it can be seen 
that the earthward flow speed is relatively constant across 
the current sheet position but the thermal pressure drops 
precipitously at that position due to the drop in ion density 
at the current sheet. The dynamic pressure behaves simi- 
larly. Thus, as the current sheet passes a fixed spatial posi- 
tion the ion momentum density drops suddenly, leading to 
the negative peak in Figure 12a. 

[ 42 ] Panel (b) of Figure 12 shows some details of the 
contributions of the individual terms on the right side of (13) 
in the pileup region. Once again, the black curve shows the 
result of a direct evaluation of the time derivative on the left 
side of (13). It can be seen that the total force density leads 
to an approximately uniform growth rate for the ion momen- 
tum in the pileup region, i.e., an earthward force density. 
The green curve shows that the contribution of the electro- 
magnetic force density is relatively minor, in agreement 
with the results of Li et al. [2011]. The red curve shows 
the contribution of the thermal pressure gradient force. 
The peak in this contribution at the leading edge of the 
pileup region reflects the energization of the reflected ions, 
leading to an increase in the temperature moment. The 
dynamic pressure gradient force dominates the remainder 
of the earthward acceleration in the pileup region. 

[ 43 ] It is important to note that (13) and the discussion of the 
preceding paragraphs relate to force balance and momentum 
change at a fixed spatial position. Li et al. [2011], however, 
studied these issues in a reference frame moving with the local 
plasma through their choice of the convective derivative 
format of the MHD momentum equation. (13) can be 
transformed to that format to obtain 

(8V ix I w SVl A re , tv ml dP ™ 

min \-dT +Vi * "&/ =*[*X+(V ( XB)J (14) 

[ 44 ] From this perspective, the dynamic pressure gradient 
term vanishes leaving the electromagnetic and thermal pres- 
sure gradient force densities, both of which are unchanged 
by the transformation. Thus, ignoring the blue curve in panel 
(b) of Figure 12, it can be seen that the thermal pressure gra- 
dient dominates at the leading edge of the pileup region, the 
electromagnetic force density peaks up at a lower value just 
in back of the thermal peak, and deeper in the pileup region 
both the electromagnetic and thermal pressure gradient force 
densities are comparable and relatively weaker. Much of the 
earthward acceleration is provided by the thennal pressure 
gradient and electromagnetic force densities at the leading 
edge of the pileup region as plasma is gathered up there by 
the oncoming disturbance. In the preceding section, we have 
seen that the ion behavior in this region is highly dynamic 
and greatly influenced by the electric field as well as the mag- 
netic field. It is not possible to reduce the electromagnetic 
term of (13) to MHD concepts generally but it is reasonable 
to think of a magnetic pressure gradient contribution at the 
leading edge of the pileup region where the field strength 
grows rapidly toward the interior of the region. To this degree, 
the DF-1 simulation results are consistent with those of Li 
et al. with an earthward total pressure gradient force density 
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Time 


Figure 19. At z=z* in Figure 18, evolution of B x (red) and 
B v (blue) with time. 

accelerating the pileup region in the simulation and with a 
reduction in the tailward total pressure gradient force density 
providing the acceleration in the magnetotail. More generally, 
without the assumption of a magnetic pressure gradient 
contribution to the electromagnetic term, it can be said that 
the DF-1 simulation results provide a somewhat generalized 
but compatible interpretation of the Li et al. results. 

6. Plasma Wave Activity in the Pileup Region 

[ 45 ] We are unaware of any reports of observed plasma 
wave activity in the upstream portion of a pileup region in 
the magnetotail. However, wave activity in this region is a 
persistent feature in our various DF simulations, including 
the DF-1 simulation. In this section, we describe the proper- 
ties of these waves in some detail and discuss the possibility 
that the wave activity is present in the observational data but 
not noticed. It appears that the wave activity might be gradi- 
ent driven by the rising pressure, field strength, and/or 
density at the leading edge of the pileup region and then 
strongly damped in back of the leading edge. There is also 


a possible relationship with the highly nongyrotropic ion 
velocity distributions in that region associated with the 
gyrophase bunched reflected ions. We have not identified 
the wave mode(s) involved. 

6.1. Current Filament Waves 

[ 46 ] One aspect of the wave activity at the leading edge 
of the pileup region in the DF-1 simulation is illustrated 
in Figure 13. The dawn-dusk current density has been 
averaged in z and then the average removed to produce the 
perturbations in the current distribution exhibited in the 
figure. The DF current sheet at x — 19 is essentially free of 
z-dependence and, so, it does not appear. The activity in the 
vicinity of the vertical dashed line, however, is unaffected 
since its z-average was originally negligible. A comparison 
with the color bar shows that the uniform blue color distrib- 
uted over the majority of the figure corresponds to J y = 0. 
Regions with yellow color contain positive current density 
and regions with darker blue or even red contain negative 
current density. The disturbance consists of filaments of 
current whose directions oscillate with varying z; the wave- 
length in this direction is approximately 2 ion inertial 
lengths. The filaments grow in front of the dashed white 
line, which is moving earthward, peak at its position, and 
decay in back of it. Each time a filament grows, it is in the 
opposite direction from the one preceding it at its position 
in z. In that sense, the pattern oscillates in the x direction 
also but the pattern is not visible at an instant in time due 
to its limited distribution in x. 

[ 47 ] We have averaged the magnitude of the current fila- 
ment waves overz. The result, at t= 10, is shown by the red 
curve in Figure 14. The blue curve shows the rise in B- in 
the pileup region, leading to the sudden jump upward at the 
position of the current sheet. The extent of the pileup region 
at this time is shown by the extent of this region of increased 
B- in front of the current sheet. It can be seen that the strength 
of the wave activity exhibits a peak that is coincident with the 
B- gradient at the leading edge of the pileup region. This is a 
position of increasing density, pressure, and flow speed 
(see Figure 4) as well. It appears possible that the wave activity 
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Figure 20. (a) Current filament wave activity in the E z electric field component at the peak wave 
amplitude position, (b) Nonideal z component, E z + V x B v . at the same position. 
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Figure 21. The passage of a current filament wave disturbance over a simulated spacecraft at three nearby 
positions. Left column: B v at (top) z = 1 , (middle) z = 0, and (bottom) z = — 1 . Right column: B x at the same 
positions. Time measured in units of co” 1 . 


is gradient driven. It should be noted, however, that the ion 
velocity distribution is significantly nongyrotropic at the leading 
edge of the pileup region due to the presence of the gyrophase 
bunched reflected ions, as shown in panels (a) and (b) of 
Figure 8. The reflected ions may be playing a role as well in 
exciting the wave activity. 

[48] The black dashed line drawn on Figure 14 is at the 

position of the wave activity peak at this time. We have 
followed the position of this peak over the course of the sim- 
ulation and have recorded the distribution in z of various 
properties of the simulated plasma along this vertical line. 
Figure 15 shows the evolution with time of the dawn-dusk 
current distribution along this line over the time interval 
t= 10 — 14.26. It can be seen that, at any fixed z, the current 
oscillates regularly for a total of approximately three wave 
periods over the time interval shown; the wave period is 
somewhat larger than 1 . At neighboring z-values, sepa- 

rated by approximately two ion inertial lengths, the waves 
are 180° out of phase; the wavelength in the z direction is 
roughly two ion inertial lengths. To evaluate these properties 
more precisely, we have carried out a 2D power spectral 
analysis of the distribution shown in Figure 15. 

[ 49 ] We have found that the 2D current filament power 
distribution contains a single peak at period 1.4 («“■') and 
wavelength 2.4 ion inertial lengths for k z > 0 and period 1.4 


(co c /) and wavelength 2.0 ion inertial lengths for K < 0. 
Figure 16 shows the power distributions on cuts through 
the peak power position for L>0. The distributions for 
k z < 0 are similar. The dimensional frequency corresponding 



Figure 22. Magnetic field magnitude measured at the 
middle row spacecraft position of Figure 21 (z = 0). Field 
magnitude including wave perturbations (red). Field magnitude 
in the absence of the wave perturbations (blue). 
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to period 1.4 is approximately 4.5 f ci , in which f ci is the ion 
cyclotron frequency. 

[ 50 ] We have calculated the electron and ion betas, averaged 
over z, for various times throughout the DF-1 simulation. 
Throughout the pileup region fi ie » 1 for all t> 7, when the 
current sheet acceleration has effectively ended. The current 
filament waves reside in high beta plasma. 

6.2. Magnetic Field Wave Activity 

[ 51 ] Figure 17 shows the evolution with time of the 
(panel a )B X and (panel b)B y field components along the 
line of maximum wave activity, constructed as discussed 
above. There is an analogous B : distribution but the distur- 
bance amplitude is relatively small; it will be neglected in 
the following. The conclusions of this section are unaffected 
by this simplification. 

[ 52 ] The field components evolve in a manner similar to 
that of the current filaments. At any z, each of the components 
oscillates in time at the same frequency as that of the current 
filaments. At neighboring z values approximately 2 ion iner- 
tial lengths apart, the oscillations are 1 80° out of phase. If this 
is the case, then there should be layers of no activity between 
the out-of-phase layers. Further, it appears as though one of 
the images is almost a copy of the other with a small shift 
in time between the two. These properties are examined fur- 
ther in the following. 

[ 53 ] Figure 18 shows the result of taking the time average 
of the magnitudes of each of the field components shown in 
Figure 17 and plotting them together as functions of z for 
comparison. Indeed, it can be seen that the activity in each 
component is layered with apparent nulls between the layers. 
We conclude from this observation that the waves in the 
magnetic field do not propagate along the background field, 
in the z direction. If the waves did propagate in the z direc- 
tion, then the nulls would be smoothed over by the time aver- 
aging done to produce the figure. 

[54] From Figure 18, it can be seen that the activity peaks 
of the two magnetic field components coincide quite pre- 
cisely. In Figure 19, we examine the relationship between 
the two components at a position z = z that coincides with 
one of the peaks; the results are typical of the behavior in 
any peak. It can be seen that the two components evolve with 
a well-defined relative phase relationship. We conclude that 
the waves are right handed elliptically (almost circularly) 
polarized with adjacent activity layers in z maintaining a 
180° out-of-phase relationship. 

6.3. Electric Field Wave Activity 

[ 55 ] The pileup region plasma just in front of the current 
sheet moves at the speed of the current sheet. The position 
of the leading edge of the pileup region moves faster as 
new, initially stationary, plasma is added. Within the leading 
edge, the plasma flow speed exhibits a steep transition from 
stationary to the current sheet speed. It is in this transition 
region that the current filament waves peak. 

[56] Panel (a) of Figure 20 shows that the E z component 
of the electric field participates in the current filament wave 
activity. The frequency, wavelength, and wave pattern are 
almost identical to those of the magnetic field components. 
Panel (b) shows the z component of the nonideal electric 
field, E : + V x B y in which both V x and B v are measured at 


the wave activity peak position; the possible V y B x contribu- 
tion is negligible. The lack of a residual signature shows 
that the current filament waves are stationary in the local 
plasma; the waves move earthward with the local plasma 
propagation speed. We have found that the flow speed 
exhibits wave activity in this region also (not shown). 
These flow speed perturbations are included in the local 
plasma propagation speed. 

6.4. Detection 

[ 57 ] The current filament waves are a persistent feature of 
our driven DF simulations but there seems to be no observa- 
tional counterpart in the magnetotail. To investigate this 
dichotomy, we simulate spacecraft observations of these 
waves by recording their signatures at various fixed positions 
as the wave front passes over. Figure 21 shows examples of 
the perpendicular B x and B y field components seen at three 
nearby positions spaced by one ion inertial length in the 
z direction. This particular group of examples has been 
chosen because the wave signatures are atypically large. 

[58] The combination of fast earthward propagation speed, 
narrow wavefield disturbance, and oscillation period commen- 
surate with the disturbance passage time leads to a short-lived 
unpredictable disturbance at any of the observation points. 
Very small spatial displacements result in the appearance of 
unrelated signatures. No clear relationship between the evolu- 
tion of the B x and B y field components is discernible. 

[ 59 ] The amplitude of these disturbances is measured in 
units of the field strength in the current sheet at the beginning 
of the simulation, whose relationship to measured field values 
in the magnetotail is remote. Notice, however, from Figure 22 
it can be seen that, in these same units, the strength of the field 
magnitude in the pileup region is close to one. Thus, in effect, 
the amplitudes of the disturbances are measured in units of the 
strength of the field magnitude in the pileup region. 

[ 60 ] Using one unit of time as a rough estimate of the 
duration of the disturbances shown in Figure 21, if the field 
strength in the pileup region were 1 nT, then the dimensional 
duration would be 10 s and if the field strength were 10 nT, 
the duration would be Is. For either case, however, the 
wave period would be commensurate with the passage of 
the disturbance. 

[61] Figure 22 shows the magnetic field magnitude 
measured at the middle row spacecraft position of Figure 21 
(z = 0). It can be seen that the entire pileup region, stretching 
from its initial growth at t — 0.4 to the sudden rise signaling 
the approach of the current sheet, is observed. The red curve 
shows the magnitude including the contributions of the current 
filament wave activity and the blue curve shows the magnitude 
neglecting the wave activity. The differences in the two curves 
lie at the leading edge of the pileup region, as expected, while 
remaining almost indiscernible. 

[ 62 ] Overall, we feel that the wave signatures shown here, 
if they were to exist in magnetotail field data, might not be 
easily recognized as indicating wave activity. 

7. Summary 

[63] We have discussed results from an ongoing driven, 
open-boundary PIC simulation study of the dipolarization 
front structure and propagation in Earth’s plasma sheet. 
Initial results from the highly idealized DF-1 model are 
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encouraging, indicating that further effort in this direction is 
warranted. Improvements, to better model the undisturbed 
plasma sheet in front of the DF, are apparent and in progress. 
These will be reported in a future publication. 

[64] Several defining properties of observed DFs are 
reproduced by the DF-1 model. Under sufficiently strong 
driving, the current sheet achieves an earthward propagation 
speed of approximately the local Alfven speed while 
retaining a coherent, numerically stable structure. The B- 
field strength in back of the current sheet and the current 
sheet intensity and thickness all self adjust under the influ- 
ence of the driver at the inflow boundary such that the current 
sheet thickness, thermal ion gyroradius at the current sheet, 
and ion inertial length evolve to approximately equal scales. 
The adjustment of the current sheet thickness to the ion 
gyroradius scale leads to an electron-ion charge separation 
in the current sheet and, consequently, to a strong normal 
electric field at, and in back of the current sheet. We have 
shown that, in the frame of the moving current sheet, this 
electric field component can be understood as due to a large 
Hall contribution but also to an almost equally large, and 
opposing, electron pressure gradient contribution. A pileup 
region forms in front of the propagating current sheet in 
which the plasma density, pressure, temperature (not shown), 
flow speed, and B : field strength rise from undisturbed 
values found far in front of the current sheet to their respec- 
tive disturbed values at the current sheet. A dawnward 
drifting region of both ions and electrons forms within 
the pileup region. Surprisingly, although the electrons 
and ions drift together they do not drift at the expected 
E x B speed. For the electrons, the difference has been 
shown to be due to an important thermal pressure gradient 
term in the momentum equation (7). We have found a similar 
explanation for the ion drifts (not shown). However, in the 
case of the ions, all of the terms in (7) play important roles 
and none can be neglected. 

[65] We have found that the DF-1 simulation exhibits ion 
reflection and acceleration in the pileup region. A model for 
the motion of the ions in the dominant B, field of the pileup 
region has been used to estimate the energy gain in the accel- 
eration process. Comparing the measured energy gain of the 
reflected gyrophase bunched ions in the DF-1 simulation to 
the prediction of the model has led to the conclusion that 
the ions are reflected by a normal electric field directed earth- 
ward at the leading edge of the pileup region in conjunction 
with the B- field there, not through an interaction with the 
current sheet and its associated magnetic field gradient. The 
reflection mechanism appears similar to that which produces 
the gyrophase bunched ions observed upstream of Earth’s 
bow shock. We have found that the highest energy ions in 
the undisturbed upstream plasma gain a factor of 2.5 - 3.5 
times their initial energy through the reflection process. 

[66] We have analyzed the contributions in the pileup 
region of the DF-1 simulation of the terms of the ion momen- 
tum equation (13) to understand the earthward acceleration 
of that region. We have found that the thermal pressure gra- 
dient dominates at the leading edge of the pileup region, the 
electromagnetic force density peaks up at a lower value just 
in back of the thermal peak, and deeper in the pileup region 
both the electromagnetic and thermal pressure gradient force 
densities are comparable and relatively weaker. Much of the 
earthward acceleration is provided by the thermal pressure 


gradient and electromagnetic force densities at the leading 
edge of the pileup region as plasma is gathered up there by 
the oncoming disturbance. It is not possible to reduce the 
electromagnetic term of (13) to MHD concepts generally 
but it is reasonable to think of a magnetic pressure gradient 
contribution at the leading edge of the pileup region where 
the field strength grows rapidly toward the interior of the 
region. To this degree, the DF-1 simulation results are con- 
sistent with those of Li et al. [201 1] with an earthward total 
pressure gradient force density accelerating the pileup re- 
gion in the simulation and with a reduction in the tailward 
total pressure gradient force density providing the accelera- 
tion in the magnetotail. More generally, without the as- 
sumption of a magnetic pressure gradient contribution to 
the electromagnetic term, it can be said that the DF-1 simu- 
lation results provide a somewhat generalized but compati- 
ble interpretation of the Li et al. results. 

[67] We have found that wave activity at the leading edge 
of the pileup region is a persistent feature of our DF simula- 
tions, including the DF-1 simulation. Wave activity growth 
coincides with the maximum gradient in density, pressure, 
field strength, etc., at the leading edge. This is also the region 
where the ion velocity distribution is highly nongyrotropic 
due to the presence of the reflected ions. Both the steep gra- 
dients and the nongyrotropic velocity distributions are candi- 
dates for exciting the waves. However, we have not identified 
the wave mode(s) involved. 

[68] The wave activity is characterized by oscillating filaments 
of current in the dawn-dusk direction. The wavelength in the 
north-south direction is approximately two ion inertial lengths, 
which is also comparable to the ion gyroradius in the pileup 
region. By tracking the position of maximum wave activity as 
it propagates with the leading edge of the pileup region, we have 
been able to determine that the waves oscillate at that position at 
approximately four times the ion cyclotron frequency. 

[69] The pileup region is characterized by high electron 
and ion betas. The current filament waves are excited in high 
beta plasma. 

[ 70 ] Perturbations in the electromagnetic field are associ- 
ated with the current filament waves. In the magnetic field, 
the waves are largely restricted to the perpendicular (to B-) 
components. The perpendicular magnetic field perturbations 
rotate about B : in a right-handed sense with layers separated 
by approximately one ion inertial length containing opposing 
perturbation vectors. The perturbation magnitudes are almost 
phase independent, leading to almost circularly polarized 
wave activity. The electric field component E- contains only 
the ideal frozen-in contribution, showing that the wave activity 
is stationary in the reference frame moving with the plasma in 
the earthward direction. 

[ 71 ] We have shown that the strength of the magnetic field 
perturbations associated with the current filament waves is 
typically a few tenths of the strength of the B- field in the 
pileup region. By placing an observation point in the wavefield 
as the disturbance propagates by, we have found that the 
combination of fast earthward propagation speed, narrow 
wavefield disturbance, and oscillation period commensurate 
with the disturbance passage time leads to a short-lived 
unpredictable disturbance at any given observation point. 
Very small spatial displacements result in the appearance of 
unrelated signatures. No clear relationship between the evolu- 
tions of the perpendicular field components is discernible. 
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[72] The presence of the current filament waves in plasma 
sheet DF observations may have gone unnoticed so far. 
Detection of the waves would be challenging. 

[73] Acknowledgments. We thank Don Fairfield for many helpful 
discussions. This research was supported by NASA’s MMS IDS grant 
NCC5-494 (MOST). 

[74] Philippa Browning thanks the reviewers for their assistance in eval- 
uating this paper. 

References 

Angelopoulos, V., C. F. Kennel, F. V. Coroniti, R. Pellat, M. G. Kivelson, 
R. J. Walker, C. T. Russell, W. Baumjohann, W. C. Feldman, and 
J. T. Gosling (1994), Statistical characteristics of Bursty bulk flow events, 
J. Geophys. Res. , 99( A 1 1 ), 2 1 ,257-2 1 ,280. 

Angelopoulos, V., et al. (2008), First results from the THEMIS mission, 
Space Sci. Rev., 747(1-4), 453-476. 

Buchner, J., and L. M. Zelenyi (1989), Regular and chaotic charged-particle 
motion in magnetotail-like field reversals. 1. Basic theory of trapped 
motion,/. Geophys. Res., 94( A9), 11,821-11,842. 

Dubyagin, S., V. Sergeev, S. Apatenkov, V. Angelopoulos, R. Nakamura, 
J. McFadden, D. Larson, and J. Bonnell (2010), Pressure and entropy 
changes in the flow-braking region during magnetic field dipolarization, 
J. Geophys. Res., 115, A10225, doi: 1 0.1 029/20 10JA0 15625. 

Eastman, T. E., R. R. Anderson, L. A. Frank, and G. K. Parks (1981), 
Upstream particles observed in the earth’s foreshock region, J. Geophys. 
Res., 56(NA6), 4379^4395. 

Gurgiolo, C., G. K. Parks, B. H. Mauk, C. S. Lin, K. A. Anderson, R. P. Lin, 
and H. Reme (1981), Non-E X B ordered ion-beams upstream of the 
earth’s bow shock, J. Geophys. Res., 56(NA6), 4415^1424. 

Gurgiolo, C., G. K. Parks, and B. H. Mauk (1983), Upstream gyrophase 
bunched ions - A mechanism for creation at the bow shock and the growth 
of velocity space structure through gyrophase mixing, J. Geophys. Res., 
55(NA11), 9093-9100. 

Hesse, M., K. Schindler, J. Bim, and M. Kuznetsova (1999), The diffusion 
region in collisionless magnetic reconnection, Phys. Plasmas, 6(5), 
1781-1795. 

Hwang, K. J., M. L. Goldstein, E. Lee, and J. S. Pickett (201 1), Cluster ob- 
servations of multiple dipolarization fronts, J. Geophys. Res., 116, A00i32, 
doi: 1 0. 1 029/20 1 0JA0 15742. 

Klimas, A., M. Hesse, and S. Zenitani (2008), Particle-in-cell simulation of 
collisionless reconnection with open outflow boundaries, Phys. Plasmas, 
15, 082102, doi: 1 0. 1 063/ 1 .2965826. 

Klimas, A., M. Hesse, S. Zenitani, and M. Kuznetsova (2010), Particle- 
in-cell simulation of collisionless driven reconnection with open 
boundaries, Phys. Plasmas, 17, 112,904, doi: 10.1 063/1 .351 0480. 
Klimas, A., M. Hesse, and S. Zenitani (2012), Particle-in-cell simulation of 
collisionless undriven reconnection with open boundaries, Phys. Plasmas, 
19, 042,901, doi: 10.1063/1.3699032. 

Langdon, A. B. (1992), On enforcing Gauss law in electromagnetic particle- 
in-cell codes, Comput. Phys. Commun., 70(3), 447-^150. 

Leroy, M. M., C. C. Goodrich, D. Winske, C. S. Wu, and K. Papadopoulos 
(1981), Simulation of a perpendicular bow shock, Geophys. Res. Lett., 
5(12), 1269. 

Li, S. S., V. Angelopoulos, A. Runov, X. Z. Zhou, J. McFadden, D. Larson, 
J. Bonnell, and U. Auster (2011), On the force balance around 
dipolarization fronts within bursty bulk flows, J. Geophys. Res., 116, 
A00i35, doi: 10.1 029/20 1 0JA0 15884. 

Mende, S. B., H. U. Frey, V. Angelopoulos, and Y. Nishimura (201 1), Substorm 
triggering by poleward boundary intensification and related equatorward 
propagation, J. Geophys. Res., 116, A00i31, doi : 1 0. 1 029/20 10JA0 1573 3. 
Nakamura, R., A. Retino, W. Baumjohann, M. Volwerk, N. Erkaev, 
B. Klecker, E. A. Lucek, I. Dandouras, M. Andre, and Y. Khotyaintsev 
(2009), Evolution of dipolarization in the near-Earth current sheet induced 
by Earthward rapid flux transport, Ann. Geophys., 27(4), 1743-1754. 


Nishimura, Y., L. Lyons, S. Zou, V. Angelopoulos, and S. Mende (2010), 
Substorm triggering by new plasma intmsion: THEMIS all-sky imager 
observations, / Geophys. Res., 115, A07222, doi:10.1029/2009JA015166. 

Ohtani, S. (2004), Flow bursts in the plasma sheet and auroral substorm 
onset: Observational constraints on connection between midtail and 
near-Earth substorm processes, Space Sci. Rev., 773(1-2), 77-96. 

Ohtani, S., Y. Miyashita, H. Singer, and T. Mukai (2009), Tailward flows 
with positive B-Z in the near-Earth plasma sheet, J. Geophys. Res., 114, 
A06218, doi: 10.1 029/2009 JA0 14159. 

Panov, E. V., et al. (2010a), Multiple overshoot and rebound of a bursty bulk 
flow, Geophys. Res. Lett., 37, L08103, doi: 10.1029/2009GL04 1971. 

Panov, E. V., et al. (2010b), Plasma sheet thickness during a bursty bulk flow 
reversal, J. Geophys. Res., 115, A05213, doi: 10.1 029/2009 JA0 14743. 

Pritchett, P. L. (2008), Collisionless magnetic reconnection in an asymmetric 
current sheet, J. Geophys. Res., 113, A06210, doi: 10.1 029/2007 JA0 12930. 

Pritchett, P. L., and J. Buchner (1995), Collisionless reconnection in configura- 
tion with a minimum in the equatorial magnetic -field and with magnetic 
shear, / Geophys. Res., 100(A3), 3601-3611. 

Runov, A., V. Angelopoulos, M. I. Sitnov, V. A. Sergeev, J. Bonnell, 
J. P. McFadden, D. Larson, K. H. Glassmeier, and U. Auster (2009), 
THEMIS observations of an earthward-propagating dipolarization front, 
Geophys. Res. Lett., 36, L14106, doi:10.1029/2009GL038980. 

Runov, A., V. Angelopoulos, X. Z. Zhou, X. J. Zhang, S. Li, F. Plaschke, 
and J. Bonnell (2011a), A THEMIS multicase study of dipolarization 
fronts in the magnetotail plasma sheet, J. Geophys. Res., 116, A05216, 
doi: 1 0. 1 029/20 1 0 JA0 1 63 1 6. 

Runov, A., et al. (2011b), Dipolarization fronts in the magnetotail plasma 
sheet, Planet. Space Sci., 59(7), 517-525. 

Schmid, D., M. Volwerk, R. Nakamura, W. Baumjohann, and M. Heyn 
(201 1), A statistical and event study of magnetotail dipolarization fronts, 
Ann. Geophys., 29(9), 1537-1547. 

Sergeev, V. A., V. Angelopoulos, S. Apatenkov, J. Bonnell, R. Ergun, 
R. Nakamura, J. McFadden, D. Larson, and A. Runov (2009), Kinetic 
structure of the sharp injection/dipolarization front in the flow-braking 
region, Geophys. Res. Lett., 36, L21105, doi:10.1029/2009GL040658. 

Sergeev, V. A., V. Angelopoulos, and R. Nakamura (2012), Recent advances 
in understanding substorm dynamics, Geophys. Res. Lett., 39, L05101, 
doi: 1 0. 1 029/20 12GL050859. 

Shiokawa, K., W. Baumjohann, and G. Haerendel (1997), Braking of high- 
speed flows in the near-Earth tail, Geophys. Res. Lett., 24(10), 1 179-1 182. 

Sibeck, D. G., and V. Angelopoulos (2008), THEMIS Science Objectives 
and Mission Phases, Space Sci. Rev., 747(1-4), 35-59. 

Sitnov, M. I., and M. Swisdak (2011), Onset of collisionless magnetic 
reconnection in two-dimensional current sheets and formation of dipolarization 
fronts, / Geophys. Res., 116, A12216, doi: 10.1029/201 1JA0 16920. 

Sitnov, M. I., M. Swisdak, and A. V. Divin (2009), Dipolarization fronts as a 
signature of transient reconnection in the magnetotail, J. Geophys. Res., 
114, A04202, doi: 1 0. 1 029/2008 JA0 1 3980. 

Xing, X., L. R. Lyons, V. Angelopoulos, D. Larson, C. Carlson, A. Runov, 
and U. Auster (2010), Plasma sheet pressure evolution related to 
substorms, / Geophys. Res., 115, A01212, doi:10.1029/2009JA014315. 

Zhou, M., M. Ashour-Abdalla, X. H. Deng, D. Schriver, M. El-Alaoui, and 
Y. Pang (2009), THEMIS observation of multiple dipolarization fronts 
and associated wave characteristics in the near-Earth magnetotail, 
Geophys. Res. Lett., 36, L20107, doi:10.1029/2009GL040663. 

Zhou, X. Z., V. Angelopoulos, V. A. Sergeev, and A. Runov (2010), 
Accelerated ions ahead of earthward propagating dipolarization fronts, 
J. Geophys. Res., 115, A00i03, doi: 1 0. 1 029/20 10JA0 15481. 

Zhou, X. Z., V. Angelopoulos, V. A. Sergeev, and A. Runov (201 1), On the 
nature of precursor flows upstream of advancing dipolarization fronts, 
J. Geophys. Res., 116, A03222, doi: 1 0. 1 029/20 10JA0 16 165. 

Zhou, X. Z., V. Angelopoulos, A. Runov, J. Liu, and Y. S. Ge (2012a), 
Emergence of the active magnetotail plasma sheet boundary from 
transient, localized ion acceleration, J. Geophys. Res., 117, A 102 16, 
doi: 10. 1 029/20 12JA0 18171. 

Zhou, X. Z., Y. S. Ge, V. Angelopoulos, A. Runov, J. Liang, X. Xing, 
J. Raeder, and Q. G. Zong (2012b), Dipolarization fronts and associated 
auroral activities: 2. Acceleration of ions and their subsequent behavior, 
J. Geophys. Res., 117, A10227, doi : 1 0. 1 029/201 2 JA0 17677. 


17 


